% read in an EDI file and compute MT app. rs. and phs.
% provision to calculate e-pre coherency 30.7.2002
% uses cohe.m for obtaining coherencey
% latest date 30.07.2002


%clear;


[file,path] = uigetfile('*.edi', 'Open EDI file');

 if file == 0,
 K = 2;
 return;
 end;

line1 = [path file];
fid = fopen(line1,'rt'); 

if (fid == -1)
	fprintf('\n\n\n\n    Error in opening %s\n',line1);
	fprintf('        ....file does not exists or path may not be correct\n');
 	return;
	end;

%-------- read in contents -----------------------------

flag = 1;
K = 0;
while flag == 1,
[A,count] = fscanf(fid,'%s',1);
if strcmp(A,'>=SPECTRASECT') == 1,
flag = 0;
end;
K = K + 1;
 if K > 5000,
	fprintf('Error in file\n');
	return;
 end;
end;

[siteID,count] = fscanf(fid,'%s',1);
[Nchannel,count] = fscanf(fid,'%s',1);
[Nfrq,count] = fscanf(fid,'%s',1);
[Nblocks,count] = fscanf(fid,'%s',2);

nfrq = str2num(Nfrq(7:length(Nfrq)));
nchan =str2num(Nchannel(7:length(Nchannel)));
[ChannelID,count] = fscanf(fid,'%f',[1 nchan]);

[head,count] = fscanf(fid,'%c',[2 1]);
for i = 1:nfrq,

if i>1,
[head,count] = fscanf(fid,'%c',[4 1]);
end;

[head,count] = fscanf(fid,'%c',[72 1]);
frq(i) = str2num(setstr(head(18:27)'));
df(i) = str2num(setstr(head(63:67)'));
% it reads only avgt 
% in future df = avgt*avgf
[data,count] = fscanf(fid,'%f',[nchan nchan]);
data = data';


HxHx = data(1,1);
HyHy = data(2,2);
HzHz = data(3,3);
ExEx = data(4,4);
EyEy = data(5,5);

HxHy = data(2,1) - j*data(1,2);
HyHx = conj(HxHy);
HxEx = data(4,1) - j*data(1,4);
ExHx = conj(HxEx);
HxEy = data(5,1) - j*data(1,5);
EyHx = conj(HxEy);
ExHy = data(4,2) + j*data(2,4);
HyEx = conj(ExHy);
HyEy = data(5,2) - j*data(2,5);
EyHy = conj(HyEy);


Zxx(i) = HyHy*ExHx - HyHx*ExHy;
Zxy(i) = HxHx*ExHy - HxHy*ExHx;
Zyy(i) = HxHx*EyHy - HxHy*EyHx;
Zyx(i) = HyHy*EyHx - HyHx*EyHy;
Den(i) = HxHx*HyHy - HxHy*HyHx;

Zxy(i)=(Zxy(i)/Den(i));
Zxx(i)=(Zxx(i)/Den(i));
Zyy(i)=(Zyy(i)/Den(i));
Zyx(i)=(Zyx(i)/Den(i));


Coh1(i) = cohe(Zxx(i), Zxy(i),data, 'Ex');
Coh2(i) = cohe(Zyx(i), Zyy(i),data, 'Ey');
[ZxxErr(i) ZxyErr(i)]  = BivError(Coh1(i),df(i),data,'Ex');
[ZyxErr(i) ZyyErr(i)]  = BivError(Coh2(i),df(i),data,'Ey');



%clear head label;
end;
fclose(fid);



%------------------------ apparent resitivity, phase----------


for i = 1:nfrq,

Phxy(i) = abs(atan(imag(Zxy(i))/real(Zxy(i))) * 180 / 3.1415926535897);
Rhxy(i) = 0.2 * (1/frq(i)) * (Zxy(i) * conj(Zxy(i)));
Phyx(i) = abs(atan(imag(Zyx(i))/real(Zyx(i))) * 180 / 3.1415926535897);
Rhyx(i) = 0.2 * (1/frq(i)) * (Zyx(i) * conj(Zyx(i)));

end;

PhxyErr = (ZxyErr./abs(Zxy))*(180/pi); % correct 31.07.2002
PhyxErr = (ZyxErr./abs(Zyx))*(180/pi);
RhxyErr = (ZxyErr.*abs(Zxy)./frq) * 0.4; % bit confusion. may be ok -dt 2.8.2
RhyxErr = (ZyxErr.*abs(Zyx)./frq) * 0.4;
%--------------------------------------------------------------
T = 1./frq;

subplot(2,1,1);
errorbar(log10([T' T']),log10([Rhxy' Rhyx']),log10([RhxyErr'/2 RhyxErr'/2]));
axis([-3 3 -1 4]);
ylabel('log ohm-m');
xlabel('log period (s)');
subplot(2,1,2);
errorbar(log10([T' T']),[Phxy' Phyx'],[PhxyErr'/2 PhyxErr'/2]);
axis([ -3 3 -180 180]);
xlabel('log period (s)');
ylabel('degrees');
title(file);

pause;
%close all;
